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ARIS  MANUAL 


1.  INTRODUCTION 

A generalized  computer  system  for  the  automated  (numerical)  integration  of  chemical 
reaction  rate  equation  has  been  developed.  The  system  is  intended  to  be  used  in  conjunction 
with  large  hydrocodes  therefore  the  integration  scheme  must  be  simple  and  efficient. 

In  order  that  such  a system  will  be  useful  for  this  purpose  it  has  to  fulfill  the  following 
criteria: 

1)  It  has  to  generate  a complete  (FORTRAN)  program  which  contains  all  "chemistry",  from 
the  user  input  to  the  system.  The  generated  program  must  be  comprehensible  (and 
always  correct).  It  must  also  contain  the  possibility  of  an  understandable  error  printout. 

2)  The  input  for  the  reaction  must  follow  the  physical  notation.  The  further  program  input 
must  be  easy  and  complete. 

3)  The  integration  schemes  provided  for  the  numerical  integration  must  be  simple  and  suited 
for  stiff  equations  with  special  features  for  asymptotic  behavior. 

The  main  features  of  the  ARIS  system  are: 

a)  Input  in  the  physical  form  of  rate  equation  for  particle  densities  and  temperatures 
which  are  changed  by  the  reactions. 

b)  Source  terms  can  be  included  with  a special  feature  if  they  represent  stiff  terms. 
ManuKript  wbinittad  Januaiy  9,  1979. 
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c)  FORTRAN  like  programming  for  initial  conditions,  rate  coefficients  and  returning  the 
integration  results. 

d)  There  is  a choice  of  three  second  order-integrators  which  are  suited  for  handling  stiff 
equations. 

e)  An  easy  communication  with  the  user  program  is  made  possible  by  allowing  users  com- 
mon blocks  (and  other  statements)  to  appear  in  the  FORTRAN  program  created  by  the  system. 

0 As  the  ARIS  system  is  intended  to  generate  automatically  a FORTRAN  program  for 
insertion  into  large  hydro  codes,  the  generated  FORTRAN  program  is  vectorizable  on  the  ASC 
computer  because  the  ARIS  system  uses  the  data  in  sequential  order. 

The  manual  describes  the  use  of  the  system.  It  gives  also  general  explanations  of  the  con- 
tent of  the  generated  subroutines  and  a short  description  of  the  integration  schemes  and  the 
integration  logic.  The  titles  to  the  different  sections  are  self  explanatory.  The  short  example 
given  at  the  end  will  make  a new  user  familiar  with  needed  input,  the  generated  FORTRAN 
program  and  the  output  provided  by  the  system. 
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II.  INPUT  SECTION 


(1)  General 


The  ARIS  system  will  generate  a FORTRAN  program  according  to  the  input  provided  by 
the  user.  There  are  basically  two  different  kinds  of  input;  REACTION  input  and  PROGRAM 
input.  The  latter  is  divided  into  five  different  sections  for  the  different  purposes  explained 
below. 

All  input  sections  start  with  a first  card  which  signifies  the  type  of  input.  This  card  begins 
with  a "#”  in  the  first  column.  The  next  such  card  (or  EOF)  closes  the  current  input.  The 
input  will  be  printed  out  as  given  except  for  this  first  card 

#RCT  Reaction  input.  Here  the  rate  equations  for  particle  densities  and  reaction  energies, 
including  conversion  factors  and  source  terms  have  to  be  given.  Rates  can  be  given  as 
parameters  or  expressions.  They  may  be  computed  either  in  ARSINT  or  ARSFUN. 

#SET  The  content  of  this  input  section  will  be  put  into  all  ARIS  subroutines  which  require 
user  input.  The  user  can  define  here  his  common  blocks  for  utilization  in  the  program 
which  sets  the  initial  conditions,  computes  rate  coefficients,  etc. 

#INT  Initial  conditions.  The  user  must  here  initialize  the  integration  variables.  The  minimum 
values  and  allowed  errors  (if  not  to  be  set  by  default). 

#FUN  User  input  to  compute  variable  reaction  rates,  reaction  energies,  conversion  factors  and 
source  terms. 

#RUN  Whereas  in  FUN  it  is  not  allowed  to  change  the  integration  variables,  the  variables  can 
be  recomputed  after  each  integration  step,  i.e.  renormalization  of  particle  densities. 

#END  Is  the  inverse  of  the  initialization.  The  user  must  specify  where  the  results  of  the 
integration  shall  be  stored  for  use  in  the  calling  program.  Functions  from  certain  terms 
can  also  be  computed  here.  They  may  be  only  available  here  because  the  content  of  the 
arrays  (ARIS  names)  will  be  destroyed  by  the  integration  of  another  point. 
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Any  one  of  these  inputs  may  be  missing,  but  the  sequence  u indicated  hu  to  be  kept.  If 
there  is  not  input  to  FUN  or  RUN  these  subroutines  may  be  missing  from  the  compiled  rou- 
tines. 

ARIS  does  simple  syntax  checks  on  the  reaction  input  and  program  input.  It  checks  also 
for  undefined  variables  and  parameters  in  their  order  of  appearance. 

If  there  is  any  syntax  error  in  the  reaction  input  the  program  will  stop  after  scanning  this 
input  with  an  error  Code  "2”. 

If  there  are  syntax  errors  on  the  program  input  or  undefined  initial  conditions,  parameters 
or  rate  coefficients  the  program  will  stop  after  scanning  the  program  input.  The  printout  will 
then  contain  in  addition  to  the  input  the  differential  equations.  If  the  user  is  developing  a pro- 
gram as  a first  step,  he  needs  only  the  reaction  input  and  can  look  at  the  printed  differential 
equations. 

(2)  Definitions  and  Conventions 

The  input  is  basically  free  format,  with  the  exception  that  a ”C”  in  the  first  column  indi- 
cates a comment  and  a "#’  in  the  first  column  designates  the  various  inputs. 

DEFINITION  ;a  STATEMENT  is  any  string  of  characters  which  does  not  contain  the  character 
The  character  ends  a statement. 

DEFINITION:  NAME  can  be  any  non  empty  character  string  of  not  more  than  eight  characters, 
which  must  not  contain  the  reaction  delimiters  ~ : , $ < > # ”.  They  may  contain  FOR- 
TRAN delimiters  -t-  / ( ).  Names  have  to  be  enclosed  in  < > with  the  exception  of  par- 
ticle densities  and  reaction  energies  in  the  reaction  input  (see  below). 

DEFINITION:  INTEGRATION  VARIABLES  are  all  names  which  appear  on  the  left  hand  side 
of  the  generated  differential  equations.  Variables  are  automatically  defined  as  ARRAYS  of 
dimension  NDIM.  They  represent  either  particle  densities  or  temperatures. 


4 


NRL  MEMORANDUM  REPORT  3954 


DEFINITION:  PARAMETERS  are  all  names  either  system  or  user  defined.  They  are  either 
scalars  or  arrays  depending  on  the  context  in  which  they  appear. 

DEFINITION:  COEFFICIENTS  are  all  rate  coefficients,  reaction  energies,  conversion  factors 
and  source  or  sink  terms.  They  can  be  either  arrays  or  scalars. 

DEFINITION:  An  EXPRESSION  is  any  FORTRAN  expression  which  may  contain 
<NAMES>. 

DEFINITION:  DECLARATION  statements:  NAME:  (parameter,  expression);  The  declaration 
statement  gives  some  directive  to  the  ARIS  system,  i.e.  PARAM:  <NAMES>;  used  in 
reaction  input  to  declare  an  integration  variable  as  a parameter. 

CONVENTION:  An  & as  index  in  a USER  name  defines  this  variable  or  parameter  as  an  array 
for  ARIS.  The  use  of  & in  INT  and  END  assumes  that  the  variable  has  at  least  the 
dimension  NTOT. 

(3)  Reaction  Inpnt 

3.1  Declarations: 

All  declaration  statements  have  to  proceed  all  other  statements  in  the  reaction  input. 

i 

PARAM:  <name  1,  ....  name  n>;  ^ 

I declares  the  variables  listed  as  parameters  which  means  they  will  be  not  integrated. 

f 

I There  can  be  any  number  of  PARAM  declarations.  The  names  may  be  repeated. 

f ■ 

j SUBST:  <name>  — 0|  < name  1>  + ...  + 02  <name  2>;  j 

I . f 


substitutes  the  variable  <name>  in  all  differential  equations.  The  coefficients  aj  must 
be  FORTRAN  expressions.  They  cannot  contain  < names  >. 
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vastly  dilTerent  for  the  linearized  equations,  i.e.  for  the  sum  and  the  difference  of 
<name  1 > and  <name  2>. 

3.2  Reaction  Statement: 

al an  1 — 61,  bm  :R2  Stl  — el,  ...  $tk  — ek; 

This  corresponds  to  the  following:  species  al,  ...,  an  are  going  to  the  species  61,  ...,  bm 
with  the  reaction  rate  coefficients  /?!  and  reaction  energies  el,  ...  ek.  And  the  inverse  reac- 
tion, the  species  61,  ...,  bm  going  to  al an  with  reaction  rate  coefficients  R2  and  reaction 

energies  —el —ek. 

Either  reaction  R 1 or  R 2 can  missing  (but  not  both) . Reaction  energies  are  optional. 

The  species  names  have  to  be  separated  by  a comma  in  the  reaction  rate  statement 
(corresponding  to  a "+"  in  the  usual  physical  rate  equations).  Each  delimiter  ends  a preceeding 
name.  The  names  here  do  no:  have  to  be  enclosed  in  brackets. 

3.3  Conversion  Statement 

<tT  : F, 

F is  the  conversion  factor  from  energy  to  temperature 

3.4  Source  Statement 

(tV  “ Source;  or  o- 1'  — %Sink; 

This  corresponds  to  an  additional  term  in  the  differential  equation  for  the  integration  vari- 
able < y>.  The  second  statement  corresponds  to  subtraction  and  is  treated  as  an  stiff  term 

(see  discussion  on  INTEGRATION).  The  ”%*  holds  for  the  entire  right  hand  side. 
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3.5  General  Explanations 

All  rate  coefficients,  energies,  conversion  factors  and  source  or  sink  terms  can  be  parame- 
ters or  expressions.  They  may  be  computed  in  AR$INT  or  ARSFUN  by  a user’s  program.  The 
ARIS  system  finds  out  if  rate  coefficients  are  constants  or  depend  (either  directly  or  indirectly) 
on  the  integration  variable.  They  will  accordingly  computed  in  AR$  INT  or  AR$  FUN  after  the 
user  input  to  such  routines.  The  differential  equation  which  result  from  the  above  defined 
reactions  and  conversions  factor  are: 

D(<aj>)IDT  ~ -RR\  + RR2  j~\ n 

D(<bj>)IDT  ~ + RR\  - RRl  y - 1,  ....  m 

with 

RR  \ ~ R\  <a\>  ...  <an> 

RR2  - R2  <bl  > ...  <bm> 

The  differential  equation  for  the  temperatures  are; 

D«tJ»IDT~fJ  (RRl  - RR2)  ej  j - 1 k 

and  for  the  source  terms  the  additional  terms 

D<  V>IDT  - source  resp.  D<  V>IDT  - -sink 
are  added. 

ARIS  counts  the  number  of  identical  appearances  of  species.  Thus  for  example 

0,0  : R - 02; 


will  result  in 
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D<Q>IDT  •-!  R <0>  <0> 

D<02>/DT  ~ R <0>  <0> 

The  reaction  statement: 

0*.  £ - 0,  £ : <0‘£>  $TE  - £; 
will  result  in 

D<0’>/DT-  <0‘£>  <0’>  <£> 

D<0>/DT  - -<0'E>  <0>  <£> 

D<TE>/DT  — -$TE <0‘E>  <0>  <£>  (conversion  factor) 

No  equation  for  <£>  will  appear  but  <£>  is  counted  as  integration  variable  unless 
specified  explicitly  as  a parameter  in  a parameter  declaration. 

(4)  Program  Inpat 

4.1  General 

The  delimiters  are  the  usual  FORTRAN  delimiters.  The  formulas  are  written  in  free 
form.  ARIS  finds  out  iteratively  by  looking  at  the  context  in  which  the  NAMES  appear 
whether  they  are  scalars  or  arrays.  It  does  simple  syntax  checks  and  tests  if  NAMES  are 
defined  in  the  sequential  order  of  appearance.  If  NAMES  are  arrays  the  loops  are  automatically 
set.  NAMES  are  only  defined  if  they  appear  on  the  left  hand  side  of  a equal  sign.  A definition 
via  a call  to  a subroutine  does  not  define  a name.  There  can  be  instances  where  a user  wants  to 
o\erwrite  either  the  automatic  definition  procedure  or  the  loop  setting.  There  are  two  declara- 
tion statements  which  will  allow  him  to  do  this. 
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4.2  Declarations 

ARRAY:  <name  1 ....  name  n>; 

The  names  so  declared  are  arrays  and  are  defined  for  the  ARIS  system.  There  can  be  any 
number  of  array  declarations. 

NOLOOP:  expression; 

There  will  be  no  loop  for  this  expression.  NAMES  which  appear  here  first  will  be  defined 
as  scalars.  (They  may  be  redefined  later  as  arrays  if  they  appear  on  the  left  hand  side  of  an 
equation).  The  index  & in  users  arrays  will  be  translated  to  1 and  the  index  for  ARIS  names 
will  also  be  translated  to  1.  If  a NOL(X)P  declaration  appears  in  #INT  the  correct  indices  for 
the  appearance  in  #END  will  be  put  into  this  expression.  A range  which  uses  NP$  in  ARSINT 
will  be  "I"  in  ARSEND.  This  feature  is  designed  for  calls  to  subroutines  but  can  also  be  used 
for  equivalences  in  #SET. 

SPLIT:  Splits  the  subroutine  ARSFUN  into  two  parts;  a)  user  input,  b)  the  computation 
of  the  rates,  conversion  factors,  etc.  This  is  useful  if  the  user  input  is  very  long. 
The  following  two  declarations  will  change  the  integration  logic.  There  are  useful  if  some  or  all 
variables  have  reached  an  symptotic  state. 

ASYMP:  (value):  <list>;  default:  value  ■■  1, 

Asymp  is  invoked  for  items  in  the  list.  If  no  list  is  given  the  default  is  all  variables. 
The  timestep  determination  ordinarily  is  based  solely  on  the  relative  difference  between 
interpolation  and  extrapolation.  With  "ASYMP"  it  is  based  on  the  minimum  of  this  and 
I creation/destruction  - l.|*value  (if  it  exists).  It  will  speed  up  the  integration  if  some 
variables  have  reached  an  asymptotic  value  before  others. 

TOT  AS:  (value):  <list>;  default  value  — 0.1,  default:  all  variables. 

This  will  terminate  the  integration  if  for  all  variables  in  the  list  | creation/destruction  • l| 
< error*  value  where  error  is  the  input  error  for  the  variable.  It  will  return  the 
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asymptotic  values.  The  integration  time  is  set  to  the  end  integration  time.  NOTE: 
This  will  only  be  implemented  if  all  specified  variables  have  creation  and  destruction 
terms.  If  several  ASYMP  or  TOT  AS  statements  are  given,  the  list  will  be  added  list. 
The  implemented  factor  is  the  last  value  given. 

PRINT;  item  I,.  . . item  n; 

The  numerical  values  of  the  above  list  will  be  printed  out.  The  current  value  for  the 
specified  point  will  be  printed.  The  items  can  be  either  a name,  or  a list  of  names,  or 
an  ARIS  FORTRAN  variable  i.e.,  (SEP  (n,  K),  K 1,3),  <E>  is  a valid  print  list. 
Each  print  declaration  gives  a line  of  output.  It  will  be  printed  at  the  beginning  of  the 
printout.  The  last  four  declarations  can  be  anywhere  in  the  program. 

4.3  Position  Statement 

As  ARIS  statements  are  free  format  there  is  sometimes  a need  to  bring  fixed  position 
statements  into  the  FORTRAN  program  i.e.  labels.  This  can  be  achieved  by  enclosing  the 
statement  in.., The  statement  itself  cannot  contain  any..*.  This  statement  does  not  end  with  ; 

Example 

"line  1 

line  2" 

The  FORTRAN  sutement  will  start  immediately  after  the  first  * the  following  line  begin- 
nings are  kept. 

4.4  Program  Inpat  deacrtption 

#SET 

The  content  of  this  section  is  put  into  each  ARIS  subroutine  to  which  is  user  input.  It  is 

placed  immediately  after  the  ARIS  generated  common.  This  feature  ia  apeciflcally  designed  to 

10 


NRL  MEMORANDUM  REPORT  3954 


I 


i 

f 


{ 


! 

allow  user  common  block  definitions  <NAMES>  cannot  appear  here,  with  the  exception  in  an 
ARRAY  declaration.  An  example  is: 

" COMMON/USER/CHEVAR  (10000,4) 

#INT 

The  user  must  here  define  the  initial  conditions,  gives  the  minimum  values  and  error  cri- 
teria. The  initial  values  can  i.e.  be  defined  by  a statement  like: 

i 

<0>  - CHEVAR  (&,  3);  I 

The  statement  will  result  in 


DO  (label)  N - 1,  Np$ 

(Label)  $V  (N,  J,  1)  - CHEVAR  (N,  3) 

where  J is  the  number  of  the  FORTRAN  variable  corresponding  to  the  name  <0>  in  the  sub- 
routine ARSINT  and  in  the  correct  statements  in  AR$END  for  the  subsequent  points. 

The  minimum  values  as  well  as  the  maximal  allowable  error  for  the  relative  difference 
between  interpolation  and  extrapolation  which  governs  the  timestep  can  be  given  for  a list  of 
integration  variables.  They  can  be  functions  of  the  integration  variables  and  parameters  but  the 
lists  cannot  be  changed.  ERRORS  and  MIMIMA  may  not  appear  on  right  hand  sides.  If  a 
name  in  the  list  is  not  an  integration  variable  a warning  will  be  given. 

If  Q represents  either  M or  E the  minima  and  errors  are  go  be  given  in  the  form: 

Q <namel  ...,  name2>  — expression; 
or 


Q<:FORALL>  — expression; 


II 
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If  the  lists  overlap  the  results  are  unpredictable. 

) The  defaults  are: 

f 

E<;FORALL>  - l.E-4; 

M<VAR>  - l.E-10  <VAR>  ; 

<VAR>  is  the  initial  value. 

The  default  holds  for  all  ERRORS  and  MINIMA  not  defined. 

#FUN 

Expressions  for  the  computations  for  the  variable  rate  coefficients,  etc.  have  to  given  here. 
If  i.e.  the  following  reaction  is  given 

N.  0+  : <NO+>  •<£>  - N+,  0 :<yV+0>  •<£>; 

Assuming  that  the  parameters  <NO+  > and  <N+0>  are  functions  of  <TE>  the  state- 
ment: 

<NO+>  - expression  (<7’£>); 

would  give  this  coefficient.  If  the  other  one  is  defined  by  a table  in  the  subroutine  TABLE  the 
’ following  two  statements  would  define  it: 

ARRAY:  <N-l-0>  ; NOLOOP:CALL  TABLE  «N+0>,.<TE>,  Np$); 

#RUN 

This  section  may  be  used  to  change  the  variables  after  one  integration-step  or  it  can  be 
used  to  rjmpute  other  functions  of  the  intergration  variables  and  parameters  used  elsewhere. 
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#END 

The  proRram  which  substitutes  back  the  integration  variaMea  or  functiona  of  theae.  It  ia 
the  inverse  to  #INT.  The  statement 

CHEVAR(ft.3)  - <0>; 


will  store  the  integrated  value  of  <0>  back  in  the  users  array. 


KLAUS  HAIN 


I 


\ 


111.  OUTPUT  SECTION 

Th»re  are  three  different  outputs  generated  either  directly  by  the  ARIS  system  or  during  a 
run.  The  last  one  consists  of  the  numerical  values  of  the  variables  and  the  rates.  This  printout 
is  according  to  the  users  needs.  It  may  be  completely  suppressed. 

(1)  Print  output 

1.1  Differential  Equations 

As  the  first  output  the  actual  differential  equations  as  generated  by  the  input  are  printed,  i.e., 
the  reaction  statement 

O.  H:  \.E-n  - HO.  1. £-24 'exp  (-/t/T); 

will  result  in  the  following 

IN1TIAL:<0>  ■=  CV  (&,1) 


d <0>/dt  - -l.E-12*<0>*<H> 

RCT  - I 

+ I.E  - 24  * exp  (-A/T)  x • <HO> 

RCT  - I 

INITIAL: <H>  - I.E  20 

d <H>/dt  - I.E  - 12  <0>*<H> 

RCT  - 1 

+ I.E-24*exp  (-A/T)*<HO> 

RCT  - I 

INITIAL:  < HO  > " undefined 

d <HO>  dt  - l.E-12*<0>*<H> 

RCT  - 1 

-l.E  - 24  * exp  (-A/T) 

RCT  - 1 

The  initial  values  for  <0>  and  <H>  have  been  defined  in  #INT,  whereas  <HO>  is 
undefined.  ARIS  will  generate  the  next  output  but  will  not  generate  a FORTRAN  program. 

1.2  Name  lists 

As  the  next  output  a table  is  generated  which  gives  the  relations  between  FORTRAN 

14 


I 


> 


NRL  MEMORANDUM  REPORT  3954 


variables  and  ARIS  names.  A general  list  is  given  in  the  next  table. 


INTEGRATION  VARIABLES  : 
PARAMETERS 
RATE  COEFFICIENTS 
SOURCE  TERMS,  etc. 

RATES 

PRODUCTION 

DESTRUCTION 

MINIMA 

ACCURACY 

REL.  ERROR  BETWEEN 

EXTRA  A INTERPOLATION  : 

INDIVIDUAL  TIME 

INTEGRATION  STEP 

OLD 

CURRENT 
INDICES  LIST 


$V(&,I,J) J - 1,2,3 

$P(&,I),  SPSCd) 
$R(S,I),$RSC(1) 

$C(&,I),  SCSC(I) 

$R(&,I)  SRSCd) 

$CR(&,I,J),  $CRSCd,J)  J -1.2 
$DS(&,I,J)— - J - 1,2 
SM(&,1),  SMSCd) 

$ER(&,I),  SERSCd) 

$EP(*,I,J)  J - 1,2 
$T(4) 

$DO(&) 

$DT(&) 

LI$(&) 


\ 

I 

I 


I 

I 

i i 


\ 


I 


i 


t 

i 

i 

t 


I - 1, . . . (number  of  specific  occurrences)  all  these  arrays  will  be  in  the  ARIS  common. 


J — 1 is  the  starting  value  for  the  integration  variables  and  the  production  and  destruction 
terms.  J - 2 is  the  extrapolated  value  of  those  variables.  J — 3 is  the  interpolated  value  of 
the  integration  variable.  For  error  terms  J - 1 is  the  old  error,  j - 2 is  the  current  error.  The 
values  for  J — 1 will  be  replaced  with  the  current  values  if  the  integration  moves  on  integration 
step  ahead. 


As  the  points  will  be  finishing  in  any  order  a list  LIS  is  kept  which  conuins  the  indices  of 
the  user  array. 
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In  general  all  ARIS  created  FORTRAN  variables  contain  a "$"  in  order  not  to  conflict  with 
users  names  defined  in  his  common  blocks.  There  are  a few  exceptions,  the  main  ones  being 
the  variables  N,  M,  L used  as  indices  in  the  ARIS  created  FORTRAN  program.  These  names 
cannot  be  used  in  the  users  common-blocks,  as  they  will  be  destroyed. 

(2)  FORTRAN  Prognm 

2.1  Subroutine  ARJMAI  (parameters) 

The  system  will  generate  a FORTRAN  program.  An  ARIS  main  subroutine  ARSMAI  is 
generated  which  in  turn  will  call  the  other  ARIS  subroutines.  It  calls  according  to  the  user  pro- 
vided input  the  different  ARIS  subroutines  whose  functions  are  described  below.  The  parame- 
ters and  the  calling  sequence  of  ARSMAI  are  described  in  the  next  chapter  "SYSTEMS 
DESCRIPTION  AND  USE." 


2.2  Subroutine  ARJINT  I 

i. 

Handles  the  initial  conditions  and  sets  up  starting  values.  User  input  to  this  routine  will 

consist  of  initialization  values  for  species  and  temperatures,  parameter  definitions  and  the  i . 

t •- 

definition  of  minimum  values  for  the  variables  and  the  maximum  allowed  relative  errors. 

These  statements  will  be  repeated  automatically  in  ARSEND  for  the  integration  of  new  points  if  , ■ 

some  of  the  points  have  been  finished.  ; F , 


2.3  Subroutine  ARSTP 

Does  the  actual  integration.  No  user  input.  There  are  three  different  integration  schemes 
possible  described  briefly  as  a comment  at  the  beginning  of  each  routine.  For  a discussion  see 
INTEGRATION  SCHEMES.  This  routine  handles  also  the  integration  logic  which  will  be 
described  later.  As  this  program  has  to  be  vectorized,  it  handles  the  logic  by  multiplication  and 
construction  of  functions  which  are  either  0.  or  I.  This  involves  more  computations  but  the 


data  flow  is  uninterrupted  and  sequential. 


This  routine  calls  ARSFUN  for  the  compuution  of 
16 


I 


I 

f 

♦ 

• i 

NRL  MEMORANDUM  REPORT  3954  [ 

) 

I 

the  variable  coefficients  and  ARSRHS  for  computation  of  the  rates,  the  total  production  and  j 

destruction  terms. 

2.4  Sabrontinc  ARSFUN 

Computes  the  variable  rate  coefficients,  reaction  energies,  conversion  factors  and  source  { 

t 

terms.  Identical  coefficients  are  eliminated,  they  will  be  computed  only  once  (identical  in  the 
literal  appearance).  User  input  is  the  program  for  these  computations.  The  system  determines 
if  rates  are  variable.  If  not,  they  are  computed  only  once  in  ARSINT  and  again  in  ARSEND  if 
there  are  needed  for  the  initialization  of  a new  point. 

2.5  SabrMtIne  ARSRHS 


Computes  the  rates  from  the  rate  coefficients  and  particle  densities.  The  total  production 
and  destruction  terms  are  then  combined  out  of  these  rates.  The  source  terms  which  may  have 
occured  in  the  reaction  input  are  added.  A similar  procedure  is  used  for  the  terms  in  the  tem- 
perature equations.  There  is  no  user  input  to  this  routine. 

2.4  Sabreatiae  AR  AUN 

Allows  the  user  to  redefine  the  integration  variables  after  one  integration  step,  i.e. 
renomalization  of  densities.  If  no  input  this  routine  will  not  be  generated. 


2.7  Sabfaatlat  ARAND 

Handles  the  logic  for  the  finishing  of  points.  If  a point  is  finished  either  in  a normal  or 

faulty  manner  the  user  programmed  resubstitution  of  the  data  for  this  point  is  executed.  The 

user  specified  calls  to  the  ARIS  print  routine  are  checked  and  the  tables  printed  out  accordingly. 

The  next  point  is  initialized  and  put  into  the  place  of  the  finished  one.  The  initialization  will  be 

repeated  as  long  as  there  are  points  left  for  which  the  integration  has  not  been  begun.  If  there 

are  no  such  points  left  the  last  point  in  the  integration  array  will  be  put  into  the  place  of  the 
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one  just  finished  and  the  number  of  integration  systems  (NP$)  will  be  decreased  by  one.  In 
this  manner  the  integration  array  is  always  densely  filled.  If  either  all  points  are  integrated  or 
the  number  of  faulty  returns  has  exceeded  the  user  specified  number  control  is  returned  to  the 
calling  program. 

2.8  Subroutine  AR5PRI 

Prints  at  its  first  call  the  time  the  ARIS  program  has  been  created  and  the  execution  time 
and  the  values  specified  in  print  declaration  (if  any).  Three  tables  can  be  chosen  by  specifying 
the  parameter  LPRI.  They  are  described  in  increasing  order  of  the  index  of  LPRI. 

1)  Variables  (densities  and  temperatures),  their  total  production  and  net  production  (pro- 
duction - destruction)  terms.  The  values  are  the  ones  at  the  end  of  the  last  successful!  integra- 
tion step  (See  INTEGRATION). 

2)  This  table  resembles  the  input  for  reactions.  The  rates  and  rate  coefficients  are 
replaced  by  their  numerical  values.  The  reaction  energies  do  not  appear.  If  a variable  is  not 
changed  by  this  reaction  the  name  does  not  appear.  The  source  terms  appear  in  a separate  table 
which  includes  the  temperatures.  The  values  of  the  rates  and  the  rate  coefficients  are  those 
used  in  the  integration.  They  are  computed  out  of  the  extrapolated  values  (see  INTEGRA- 
TION). 

3)  All  rates  and  source  terms  will  appear  under  the  variable  name  with  the  reaction 
number  which  created  them.  As  these  rates  are  the  same  for  all  species  involved  in  a reaction 
this  will  generate  much  more  output  than  Table  2.  The  numerical  values  for  the  extrapolated 
and  interpolated  variables  are  given.  For  the  temperatures  only  the  source  terms  appear.  The 
rates  and  the  total  production  and  destruction  terms  are  those  currently  used,  employing  the 
extrapolated  values.  The  destruction  terms  may  differ  from  the  sum  of  the  destruction  terms  as 
destruction/variable  rather  than  the  destruction  term  itself  is  used  by  the  integration  routine. 
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IV.  ARIS  SYSTEM  AND  ITS  EXECUTION  j 

In  order  to  use  the  ARIS  system  the  user  must  first  invoke  the  execution  of  the  ARIS 
program,  which  is  described  below  for  the  ASC.  After  a successful  execution  and  a subsequent 
FORTRAN  compiliation  he  has  to  call  the  subroutine  ARSMAI  (parameters)  from  his  program. 

The  parameters  are  described  in  the  next  section. 

(1)  ARIS  Execution  (JSL)  j 

I 

The  following  job  control  statement  holds  only  for  the  Texas  Instruments  ASC.  It  will 
not  be  difficult  to  provide  a similar  procedure  for  any  other  computer.  In  order  to  execute  the 
"ARIS”  program  one  proceeds  in  the  following  way:  First  a MACRO  ASSIGN  Statement  and 

f 

then  an  actual  EXECUTION-Statement.  | 

/ MACASG  SDUM,USERCAT/D77/BOO/HAINKl/TERMI/MACROX  ' 

/ ARIS  (ARISDAT-INPUT,ARISOUT-FPROG,RES-PRINT,DIM-NDIM,OPT-N, 

COMMON  - N AME.CHEM  - C) 

All  parameters  are  optional  and  have  default  values. 

ARISDAT  - 
INPUT  name  of  input  file. 

DEFAULT:  input  cards  following  the  / ARIS  control  card 
ARISOUT  - 

FPROG  name  of  FORTRAN  program  Ale. 

DEFAULT:  ARISOUT 

RES  - 

PRINT  name  of  print  Ale. 
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DEFAULT:  System  print  out. 

DIM  - 

DIM— number  of  points  to  be  integrated  simultaneously.  The  storage  space  is  approxi- 
mately: Storage— NDIM  (12  variables-l-3  number  of  reactions+6  +user  defined 
arrays) 

DEFAULT:  200 

OPT  - n 

n - 1,2,3  invokes  the  different  integrators.  For  discussion  see  below. 

DEFAULT- 1 

COMMON  - NAME 

NAME  common  area  for  ARIS  scratch  storage  space. 

DEFAULT:  Blank  common 

CHEM  - C 

Changes  the  first  character  of  the  names  of  the  generated  subroutines.  All  subroutines 
will  start  with  CRS 

DEFAULT  - A 

t 

In  order  to  enable  the  user  to  determine  at  which  time  the  subroutines  have  been  gen- 
erated all  subroutines  (including  ARSMAI)  have  as  a comment  the  time  at  which  they  have 
been  created.  This  will  give  an  indication  of  the  connection  with  the  input  if  a complete  recom- 
pilation is  not  done  every  lime. 

If  the  logic,  which  means  the  reactions  as  such,  the  number  of  variables  and  parameters 
including  their  order  of  appearance,  has  not  been  changed,  only  those  routines  in  which 
changes  occur  have  to  be  recompiled. 


I 

I 
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The  user  can  also  make  changes  directly  in  the  generated  subroutines  by  using  some  edi- 
tor programs.  The  meaning  and  use  of  indices,  the  connection  between  ARIS  names  and  FOR- 
TRAN variables  has  been  explained  in  the  output  section. 

(2)  Program  Execatloa  (calliag  soqaencc) 

In  order  to  execute  the  compiled  ARIS  subroutines  the  user  has  to  call  the  main  subrou- 
tine ARSMAI  with  several  parameters  which  are  described  below. 

CALL  AR$MAI(NTOT,TIMEA,TIMEB,DTO,DTMIN,LER,NEMAX.NER.NF,LPRI)  where; 

NTOT: 

Total  number  of  points  to  be  integrated.  This  number  can  be  greater  or  less  than  the 
number  of  points  (NDIM)  to  be  integrated  simultaneously. 

TIMEA: 

Starting  time  for  integration. 

TIMEB: 

End  time  for  integration 

DTO: 

Guess  for  the  first  integration  step.  The  integration  scheme  will  adjust  the  integration 
stepsize  automatically  according  to  the  desired  accuracy.  As  the  integration  step  is 
always  to  be  assumes  positive,  different  signs  indicate  the  following: 

DT0>0 

Array  input  (dimension  NTOT).  The  last  integration  stepsize  is  returned,  such  that  for 
a following  integration  the  starting  integration  stepsize  is  correct. 

DT0<0 

Scalar  input.  No  return  values. 
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CAUTION: 

The  user  has  to  make  sure  that  his  time  increments  do  not  fall  below  the  roundoff 
error.  (DT<l.E-6  TIME  A).  If  necessary  he  has  to  set  TIMEA  — 0 and  replace  TIMEB 
by  the  difference  adding  TIMEA  where  absolute  times  are  needed. 

DTMIN: 

The  minimum  timestep  allowed.  The  ARIS  program  will  increase  NER  by  one  for 
every  point  for  which  the  integration  step  falls  below  DTMIN.  The  index  of  these 
faulty  points  will  be  put  into  LER(NER).  The  integration  will  be  discontinued  for  this 
point. 

LER: 

An  array  (dimension  NEMAX)  will  contain  the  number  index  faulty  points,  if  any. 

NEMAX: 

The  maximum  number  of  faulty  returns  allowed,  it  is  set  by  the  user. 

NER: 

The  number  of  faulty  points.  If  NER>NEMAX  the  program  will  return  to  the  calling 
program.  The  remaining  points  are  neither  finished  nor  are  their  return  values  set. 

NF: 

The  total  number  of  right  hand  side  (function)  calls.  This  will  give  an  indication  of  the 
speed  of  the  integration. 

LPRI: 

(dimension  3)  controls  the  Printout. 

LPRKJ)  - 0 : 
no  print  out 
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LPRI(J)  - N : 

for  n points  which  finish  first  either  regularly  or  faulty  the  selected  table  is  printed. 
LPRI(J)  - -N  : 

up  to  n faulty  points  are  printed  out  but  as  ARIS  returns  after  NEMAX  faulty  points. 
No  more  than  NEMAX  points  will  be  printed. 

LPRKl): 

Variables,  total  production  and  net  production  terms.  The  values  are  those  used  for  the 
extrapolation  in  the  next  integration  step.  They  are  the  old  values  if  the  integration 
step  was  not  completed  (see  INTEGRATION). 

LPRI(2): 

Gives  in  addition  to  LPRI(I),  the  rates  and  rate  coefficients  sorted  according  to  reac- 
• tions.  This  table  will  resemble  the  reaction  input,  but  replaces  the  names  with  numeri- 

cal values.  The  source  terms  are  given  in  separate  tables.  Whereas  the  rates  are  given 
I only  for  the  densities,  the  source  terms  include  the  temperatures.  The  rates  and  source 

terms  given  are  those  used  for  the  current  integration  step. 

I LPRI(3): 

Prints  the  variables,  total  production  and  destruction  terms,  rates  and  source  terms 
f sorted  according  to  species.  The  values  are  those  used  in  the  last  integration  step  com- 

puted by  using  the  extrapolated  values.  The  total  destruction  terms  may  differ  consid- 

4 

erably  from  the  sum  of  the  rates  if  an  integration  step  has  not  been  completed  and 
^ slighty  in  a normal  integration  step. 

I 

For  the  variables  the  extrapolated  and  interpolated  values  are  given.  By  comparing  these  two 
values  the  user  may  have  the  possibility  of  finding  which  variable  determines  the  timestep. 
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V.  INTEGRATION 

(1)  GcneniS  Remarks 

As  the  system  is  intended  to  be  used  in  conjunction  with  large  hydrocodes  the  integration 
formulas  should  be  simple  and  fast. 

The  numerical  integration  should  be  able  to  handle  very  stiff  equations  reasonably  well 
including  temperature  equilibration.  The  program  should  never  stop  unless  there  are  user  mis- 
takes in  the  coefficients  which  cause  overflow  or  similar  conditions. 

The  accuracy  need  not  to  be  too  high.  It  therefore  seems  appropriate  to  use  a second 
order  integrations  scheme  with  predictor  and  corrector  similar  to  CHEMEQ  used  intensively  in 
the  U.  S.  NAV.  RES.  LAB  (1). 

The  integration  logic  is  simple,  namely:  If  the  relative  error  between  extrapolation  and 
interpolation,  or  for  the  asymp  option:  the  minima  of  this  eiior  and  the  ^roduction/destruction 
- l|  is  exceeds  the  error  criterium  by  a factor  of  four  the  timestep  will  be  reduced  by  a factor  of 
four  and  the  step  will  be  repeated.  If  the  error  exceeds  the  given  error  criterium  up  to  a factor 
of  four,  the  timestep  will  be  halved  but  the  integration  is  advanced;  if  the  error  is  less  than  1/4 
the  error  criterion  the  timestep  will  be  increased  by  *yJ2.  There  are  no  iterations  because  it  is 
usually  more  profitable  to  proceed  with  a smaller  timestep  than  to  spend  time  in  iterations. 

Because  the  timestep  is  reduced  very  quickly  the  user  is  advised  to  use  as  a first  guess  a i 

step  which  is  more  likely  too  high  than  too  low.  For  repeated  integrations  the  individual  I 

* t 

timestep  as  returned  after  the  integration  should  be  used.  | 

< 

There  are  three  different  integration  schemes  provided.  They  differ  in  the  handling  of  the 
stiff  terms.  This  will  be  explained  below.  The  last  one  (OPT  <-3)  consists  of  the  formula  used  | 

by  CHEMEQ  (without  asymptotics).  Generally  speaking  all  three  schemes  give  essentially  the  j 

I 

same  results  in  about  the  same  time.  The  la«t  OPT~3  will  run  a little  bit  faster  if  the  equations  ' 
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are  not  too  stiff.  OPT  •>  1 insures  the  convergence  of  the  integration  in  all  cases  and  should  be 
used  for  very  stiff  equations. 

(2)  Integration  SehenM 

The  set  of  differential  equation  for  chemical  reactions  are  in  general  of  the  form 


(2.1) 


where  d-y  is  the  destruction  term  and  c the  production  term.  Assuming  d and  c constant  one 
can  integrate  directly  and  obtains 


yit)  - y(o)  e"‘'  '+  (1  - 0 

a 


(2.2) 


As  discussed  before,  the  ARIS  system  uses  a second  order  integration  scheme.  The  different 
options  are  distinguished  by  how  they  approximate  the  factor 

/ - e-"  ' (2.3) 

during  a timestep.  In  general,  one  can  replace  this  factor  by  any  fractional  polynomial  so  long 

as  it  represent  the  exponential  function  to  at  least  the  second  order  of  d-t.  As  the  system 

should  handle  stiff  equations  where  in  one  time  step 

d t » 1 (2.4) 

the  factor  should  go  to  zero  for  large  t.  For  second  order  accuracy  of  /can  be  represented  by 


/- 


1 — a(d't) 


(2.5) 


1 + (1  - a)  d t+fi(d  t)^ 

Then  1// should  correspond  to  the  exponential  function  to  at  least  the  second  power  in  d-i.  It 


follows  that 


(2.6) 


The  different  option  are  characterized  by 

opf  * I a -•  0 ^ “ 1/2 

opt -2  a - 1/3  P - 1/6 
opt  — i 0-1/2  /5  - 0 

opt  - 2 represents  the  exponential  function  exact  to  the  third  power,  opt  ~ 3 corresponds  to 


(2.7) 


CHEMEQ. 
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The  integration  performs  first  an  extrapolation 

>>^(1  - a{d„  t))  + r c„(l  + P d],-  t) 


ye 


Then  an  interpolation  given  by 


1 + (1  - a)  • / + /3(</„  • /)2 


y.  - 


1 -y  ■ a ■ {d„  + d,) 


+ y • (fo  + + P 


1 + y • (1  - a)  • (d,  + d,)+pU  d,y 


which  will  give  for  /3  ^ 0 for  stiff  equation 


y,(asymp)  - 


(2.8) 


(2.9) 


(2.10) 

The  formula  used  in  the  different  options  are  given  as  a comment  in  the  corresponding  integra- 
tion r‘’'Jtine  (ARSSTP). 


If  the  declaration  PAIR:  <namel,  name2>;  was  used  the  transformation  to  the 
corresponding  eigenvectors  are  performed  before  the  integrationstep  backtransformation  after 
the  extrapolation  and  the  intrapolation.  The  matrix  elements  are  computed  only  once  during  at 
timestep. 


(3)  Integratlen  Legk 


The  relative  error  dfiy)  which  determines  the  timestep  is  defined  as  follows: 

d/iy)  - \y,  - y\ly- 

or 

df(y)  - min(|y,  - y\lyQ\  factor  x|*  x ye/c^  - l|) 
if  ASYMP  has  been  specified  (and  c,  exists).  The  prescription  for  determining  the  timestep  is: 


rf/(x)>4*error  : 

<ft-yrff 

no  advance 

4*error 

> 

dff(x)  > error  : 

advance 

error 

> 

d/(x)  > yerror  : 

dt^dt 

advance 

yerror 

> 

dfiy) 

1 

1 

advance 

I 
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Three  FORTRAN  arrays  (SFR)  are  constructed  which  are  either  0 or  1 depending  in  what 
range  of  dfiy)  lies.  The  advanced  and  the  old  variable,  prediction  and  terms  are  multiplied  by 
a combination  of  these  $F|l  arrays  such  that  the  above  procedure  is  applied.  The  same  tech- 
nique is  used  to  determine  the  new  timestep. 

The  starting  of  one  point  is  achieved  by  setting  the  integration  step  length  to  zero  which 
will  automatically  increase  the  timestep  by  the  factor  4>^.  Therefore  this  factor  is  applied  to 
the  user  given  guess.  Care  is  also  taken  to  reach  the  exact  TIMER. 

If  TOTAS  has  been  specifled  then  if  for  all  \d.ylc  — l|  <error  factor  the  time  is  set  to 
TIMER  and  so  the  integration  will  be  terminated. 
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VI  EXAMPLE 
1)  INPl/r 

This  example  Is  a small  portion  of  an  ionospheric 
problem.  It  shows  the  basic  features  of  the  APIS  system. 
All  neutral  species  are  declared  as  parameters  because  one 
knows  beforehand  that  they  would  not  change  appreciably  in 
this  altitude  region.  The  common  defined  in  set  are  for  the 
whole  problem  so  they  are  more  complicated  as  needed  for  this 
example.  The  function  LOOKTE  is  a table  lookup  for  all 
FORTRAN  variables  in  the  common  block  RTLOOK. 

The  input  is  listed  here  in  a form  as  returned  by  the 
ARTS  system.  Ihe  only  difference  between  this  table  and  the 
data  cards  is  the  appearance  of  the  sequence  numbers  and  the 
starting  card  ( indicated  In  the  comment). 

APIS  INPOT  ••  REACTIONS  **  03/03/77  10:55 

C #RCT 

C 

C SHORT  EXAMPLE 

C 

RCT  = 0 PARAM:<0,02,01D>;  PARAM:<N2D>; 

C DECURE  ALL  NEUTRALS  AS  PARAMETERS 

C 

RCT  s 1 N+  ,02  : <N4O2>»1.30  = N2D,02+ 

$TI=.045 

RCT  = 2 N+  ,02  ; <N402>».6^  s 0 ,N0+ 

ITIs3.32; 

RCT  = 3 N0+  : ADR(&,1)«<E>  e N2D  ,0 

: 1.il2E-12»EXP(-.i<»<TI»; 

RCT  = 4 02+  :ADR(4,2)«<E>  a 0 ,01D 

$TEs1 .5»«TE>-<TI»+1 .25; 

C 

C CONVERSION  FACTOR  FOR  TI  AND  TE 

C 

RCT  s 5 OTI  : 2./(3.«<E»;  WE  : 2./(3.»<E»; 

C 

C EQUILIBRATION  BETWEEN  TE  AND  TI 

C 

RCT  *7  WI  = TEQ0NE(&)»<TE>;  WE  s TEQ0NE(4)«<TI>; 

RCT  s 9 WE  =JTEQ0NE(4)»<TE>;  WI  s»TEQ0NE(4)«<TI>; 

C 

C END  OF  REACTION 

C 
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APIS  INPUT  •»  SET  »»  03/03/77  10:55 


SET  = 
SET  = 
SET  = 

SET  = 

SET  s 
SET  s 
SET  = 
SET  = 


COMMON  FOR  OUTSIDE  VARIABLES. 

COMMON/ATMOS/  C0(151 ) ,CN2(151 ) ,CN(151 ) 
. ,TN(151),VMID(151),CNP(151),COP(151) 

. ,TMPI(151),TMPE(151)" 

COMMON  FOR  MOLECULES(IONISED) 

COMMON /MCLI0N/C02P  (151)  ,CN0P  (151)'' 
COMMON  FOR  LOOKUP  TABLE 

COMMON/RTLOOK/  AN(100,2) ,AN2D(100,2) 

. ,A0(100,2),A01D(100,2),AN2(100,1,) 

. ,ANP(100,1),AOP(100,1),ADR(100,2) 

. , AEX(100,12),TEQONE(100)" 

END  OF  COMMON 


ARIS  INPUT  ••  INITIAL  »»  03/03/77  10:55 


#INT 

INITIALIZE  VARIABLES  AND  NEUTRAL  DENSITIES. 


= 1 <0>  = 00(4);  <02>  = 002(4); 

s 3 <01D>  = 0.  ; <N2D>  = 0.  ; 

s 5 <N+>  = CNP(4);  <N0+>  = CN0P(4); 

s 7 <02+>  s C02P(4); 

0 INITIALIZE  TEMPERATURES 
= 9 <TI>  s TMPI(4);  <TE>  s TMPE(4); 

0 SET  MINIMUM  AND  ERRORS 
0 

0 MINIMUM  TEMPERATURES  ARE  SET  BY  DEFAULT 
0 

= 11  M<H>  = 1.E-6»(<0>+<02>+<N2D>+<01D»; 

= 12  M<E>  = <02+>+<N0+>+<N+>; 

= 13  M<02+>  = AMAX1«02+>,<H»; 

= 14  M<N0+>  = AMAX1«N0+>,<H>); 
s 15  M<N-».,N0+,02+>  = AMIN1<1.E-6»<E>,<H»; 

= 16  E<:FORALL>  s 1.E-6; 

0 END  OF  INITIAL  CONDITONS 
0 


ARIS  INPUT  ••  FUNCTIONS  ••  03/03/77  10:55 


IFUN 

LOOKTE  AS  A TABLE  LOOKUP  AS  FUNCTION  OF  TE 
FOR  ADR  AND  TEQONE 
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FUN  = 1 NOLOOP:  CALL  LOOKTE«TE>,NP$); 

C COMPUTE  THE  OTHER  RATE  COEFFICEffTS 
C AND  <E>  FOR  THE  CONVERSION  FACTOR 
FUN  z 2 <E>  z <N+>  + <02+>  + <NO+>; 

FUN  z 3 <LTI>  z ALOG«rE»; 

FUN  z 4 <H>  = SIGN(.5,<TI>-.396); 

C COMPUTE  SIGNUM  FUNCTION  FOR  COMPUTATION 
C OF  RATE. 

C 

FUN  z 5 <N402>  =2.8E-10»(.5-<H» 

FUN  z 5 4.8E-12»(.5+<H»»EXP(.57»<Ln»; 

0 END  OF  FUNCTION  EVALTIONS 
C 

ARIS  INPUT  ••  END  ••  03/03/77/10:55 


C lEND 

C RETURN  INTEGRATED  VALUES 
C 

END  z 1 CNP(4)  = <N+>;  C02P(&)  z <02+>; 

END  z 3 CNOP(&)  = <N0+>; 

C END  OF  INPOT 
C 

2)0in‘PUT 

The  ARIS  program  will  print  the  error  message  (if  any) 
after  the  line  in  which  the  error  has  occurred.  There  are  a 
few  exceptions.  The  main  exception  is  that  the  coefficients 
will  be  checked  for  FORTRAN  syntax  errors  after  establihing 
the  differential  equations.  Therefore  any  error  (warning) 
messages  will  appear  immediately  after  the  differential 
equations. 

ARIS  OUTPUT  DFFERENTIAL  EQUATIOKS  03/03/77  10:55 

•••  RATE  EQUATIONS  ••• 

EQU  z 1 INITIAL  <N-f>sCNP(&) 

D(N+  )/DT  z -<N>»<O2>»«N^2>»1.30)  RCT  z 1 

-<N+>«<02>»«N402>».64)  RCT  z 2 

M»e 

EQU  z 2 INITIAL  <024.>zC02P(&) 

D(02^  )/DT  z ♦<H+>«<02>»«N-*02>»1.30)  RCT  z 1 

-<02^>»(ADR(A,2)»<E>  RCT  z 4 

•••• 

EQU  z 3 INITIAL  <N(>4->zCN0P(&) 
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D(NO+  )/DT  = +<N+>«<02>»(<N402>».6i|)  RCT  = 2 ^ 
+<N 2D>»<0>»n . 42E-1 2»EXP (- . M»<TI » 

-<N0+>«(ADR(4,1)»<E»  RCT  = 3 


«««• 


TEMPERATURE  EQUATIONS  *»» 

EQU  = 4 INITIAL  <ri>  = TMPI(4) 

D(TI  )/DT  = (2./(3.»<E»)«(  FAC  = 0 

♦<N+>«<02>»(<N-»02>»1 . 30) 

•(.045) 

+<N+>»<02>«( <N402>* . 64 ) 

•(3.32)) 

+(TEQ0NE(4)^<TE» 

J(TEQ0NE(4)^<TI> 

EQU  = 5 INHIAL  <TE>  = TMPE(4) 

D(TE  )/DT  = (2./(3.^<E»)^( 

♦<02+>^(ADR(4,2)^<E»^ 
(1.5^«TE>-<ri»+1.25)) 
♦(TEQ0NE(4)^<TI» 

»(TEQ0NE(4)^<TE» 

•••••  WARNING:  USE  CF”  "A”  IN  VARIABLE  RATE  COEFFICENTS 
MAY  CAUSE  ERRORS 

It  follovB  part  of  the  namelist  in  order  to  show  the 
connection  between  ARIS  names  and  FORTRAN  variables. 

ARIS  OUTPUT  NAMELIST  03/03/77  10:55 

VARIABLES 


N-t- 

$V(4,3,1) 

N04> 

$V(4,2,1) 

02+ 

$V(4,1,1) 

TE 

$V(4,5,1) 

TI 

$V(4,4,1) 

PARAMETERS 

E 

$P(4,4) 

H 

$P(4,5) 

LTI 

$P(4,6) 

N-*02  $P(4,7) 

NO 

4PSC(1) 

N2D 

$PSC(3) 

0 

$P(4,1) 

0+ 

$P(4,3) 

01D 

4SC(2) 

02 

$P(4,2) 

CREATION  TERMS 

H* 

VOID(O) 

N04- 

$CR(4,2) 

02+ 

$CR(4,1) 

TE 

$CR(A,4) 

TI 

ICR (4, 3) 

RCT  = 1 

RCT  = 2 
EXT  = 0 
E)Cr  = 0 


FAC  = 0 

RCT  = 4 
E)Cr  = 0 

E}cr  = 0 
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DESTRUCTION  TERMS 

N+  $DS(&,1)  NO+  $DS(4,3)  02+  $DS(4,2) 

TE  $DS(4,5)  TI  $DS(4,4) 


RATE  COEFFICIENTS 


RCT 

1 

$C(4,1) 

FUN 

VOID(O) 

RCT 

2 

4C(4,2) 

FUN 

VOID(O) 

RCT 

3 

$C(4,3) 

FUN 

$C(4,4) 

RCT 

4 

$C(4,5) 

FUN 

VOID(O) 

FUN 


3)  FORTRAN  PROGRAM 

In  the  following  a part  of  the  created  FORTRAN  progran 
is  given.  All  loops  are  over  n from  1 to  NP$  are  ommitted  in 
order  to  show  more  clearly  the  created  program. 

SUBROUTINE  AR$FUN(J) 

C ++♦  ARIS  EXECUTION  HAS  BEEN  AT:  03/03/77  10:55 
Aris  common 
User  common 

CALL  L0KTE($V(1,5,J),NP$) 


table  lookup  as  a function  of  TE 


$P(N,4)=$V(N,3,J)+$P(N,3)+$V(N,1,J)+$V(N,2,J) 
$P(N,6)  =AL0G($V(N,4,J)) 

$P(N,5)  =SIGN(.5,$V(N,4,J)-.396) 

$P(N,7)  s2.8E-10»(.5-$P(N,5))+4.8E-10‘»EXP(.57«P(N 
,6))»(.5+$P(N,5)) 


The  users  program  ends  here.  It  follows  the 
computations  for  the  coefficients.  One  can  note  that  the 
conversion  factor  $C(N,7)  appears  only  once,  as  well  as  the 
source  and  sink  terms  for  the  temperatures 
$C(N,8),$C(N<9). 

$C(N,1)  5 $P(N,7)»1.30 
$C(N,2)  s $P(N,7)».64 
$C(N,3)  = ADR(N,1)»$P(N,4) 

$C(N,4)  s 1.42E-10«EXP(-.4»$V(N,4,J)) 

$C(N,5)  = ADR(N,2)»$P(N,4) 


I 


The  coefficients  for  the  energy  changes  follow  now 

32 


■Hiiam 


f 


( 

e 

t 

t 

i 


\ 

NRL  MEMORANDUM  REPORT  3954 


The  two  constants  for  the  first  two  reactions  are  set  in 
AR$INT 

$C{N,6)  = 1.5»($V(N,5,J)-$V(N,4,J))+1.25 
The  conversion  factor 


$C(N,7)  = 2./(3.*P(N,4)) 

Source  and  sink  terms 

$C(N,8)  = TEQ0NE(N)»$V(N,5,J) 
$C(N,9)  = TEQ0NE(N)»$V(N,4,J) 
RETURN 
END 


In  the  following  subroutine  AR$RHS  the  rates  and  total 
creation  and  destruction  terms  are  computed. Part  of  this 
routine  is  given 

SUBROUTINE  AR$RHS(J) 

C +++  ARIS  EXECUTION  HAS  BEEN  AT;  03/03/77  10:55 
aris  common 


$R(N,1)  = $C(N,1)»$V(N,3,J)»$P(N,2) 
$R(N,2)  = $C(N,2)»$V(N,3,J)»$P(N,2) 
$R(N,3)  = $C(N,3)*$V(N,2,J) 

$R(N,4)  = $C(N,4')»$PSC(3)»$P(N,1) 
$R(N,5)  = $C(N,5)»$V(N,1,J) 


and 


It  follows  the  computation  of  the 
destruction  terms. 


$CR(N,1,J)  = +$R(N,1) 
$CR(N,2,J)  = ♦$R(N,2)+$R(N,4) 


total  creation 


The  energy  changes  for  the  temperatures  are  regarded  as 
creation  term. (Only  TE  is  given) 

$CR(N,3,J)  = $C(N,7)»(+$R(N,1)»$CSC(1) 
■4R(N,2)»$CSC(2))+$C(N,8) 

The  destruction  terms  are  constructed  similar  to  the 


i 

\ 
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creation  terms.  They  are  divided  by  the  variable  itself  at  : 

the  end  of  the  subroutine.  | . 

$DS(N,1,J)  = $DS(N,1,J)/$V(N,3,J) 

RCTURN 

END 

In  order  to  understand  the  integration  logic  the 
essentials  of  the  integration  routine  is  given  below. 

SUBRUTINE  AR$STP(N$F) 

C +++  ARIS  EXECUTION  HAS  BEEN  AT:  03/03/77/  10:55  i 

C 

«»•  OPT  = 1 

C YE  = (Y4DT»CR»(1.+.5»DT»DS))/(1.4DT»(DS»(1.+.5»DS)) 

C 

C n = (Y+.5<^T»(CR4CRE»(1.4DT«^)SE))) 

C /(1.+.5»DT«(DS+DSE*(1.+DT«DSE)))  f 

c ! ■ 

Aris  common  . 

Factors  for  the  errors  are  computed  first 

REAL  $FAC(3)/.25,1.,J4./ 

DO  8000  K=1,3 

8000  $EF(1,K)  = $ERSC(1)»$FAC(K)  , • 

For  all  scalar  error  criteria  *■  . 

for  all  destruction  term  compute  i . 

I 

$H(N,L,1)  = $DT(N)«$DS(N,L,1) 

$H(N,L,2)  = 1.+.5«$H(N,L,1)  [ 

I 

Compute  the  extrapolated  variables.  Note  the  absence  of  I . 

the  creation  term  in  the  first  formula. (Only  two  variables  ! ’ | 

are  given). 

$V(N,3,2)  = $V(N,3,1) 

/(1.4$H(N,1,1)»$H(N,1,2))  f ; ' 

$V(N,1,2)  = ($V(N,1,1)+$DT(N)»$CR(N,1,1)»$H(N,2,2))  . f; 

/(1.4W(N,2,1)«$H(N,2,2))  I' 
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compute  the  maxima 

$V(N,L,1)  = AMAX1($V(N,L,1),$M(L)) 

CALL  AR$FUN(2) 

CALL  AR$RHS(2) 


I 

I 

I 


f 


the  interpolated  values  are  computed  next  after  the 
computation  of  the  right  hande  sides  in  AR$FUN  and  AR$rhs 

$DTH(N)  = .5»$DT(N) 

$H(N,L,1)  = 1.  +$DT(N)«$DS(N,L,2) 

$V(N,3,3)  = $V(N,1,1) 

/(l.+$DTH(N)»($DS(N,1,1)+$DS(N,1,2) 

*$H(N,1,1))) 

$V(N,1,3)  = ($V(N,1,1)+$DTH(N)« 

($CR(N,1,1)+$CR(N,1,2)»$H(N,2,1))) 

/ ( 1 . +$DTH (N )• ($DS (N , 2, 1 )«$DS (N , 2, 2) 
•$H(N,2,1))) 


• • • • 


compute  maxima 

compare  for  all  variables  extrapolation  and  intrapolation 

$EP(N,L,2)=.5»$EP(N,L,1)+ 

ABS($V(N,L,2)-$V(N,L,3))/$V(N,L,1) 

compute  the  signim  functions  for  the  relative  errors 
and  sum  them  over  all  variables. 

$FR(N,K)  = SIGN(1.00001,$EF(1,K)-$EP(N,3,2)) 

$FR(N,K)  = $FR(N,K)  + 

SIGN ( 1 . 00001 , $EF( 1 , K)-$EP(N , 1 , 2) ) 


Then  compute  the  total  sign  function. (Note  that  5.  is 
the  nunber  of  variables) 


$FR(N,K)  = .5  + SIGN(.5,$FR(N,K)-5.) 


The  following  table  shows  the  values  of  $FR  and 
their  connection  with  the  timestep  changes  and  advances. 
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K = 

1 

2 

3 

$FR 

1 

1 

1 

DT=DT»2»»(1/4) 

advance 

$FR 

0 

1 

1 

DT=DT 

advance 

$FR 

0 

0 

1 

DT=DT/2 

advance 

$FR 

0 

0 

0 

DT=DT/4  no 

andvance 

The  formula  for  the  new  advanced  timestep 
then  is  the  following 

$DO(N)  = $D0(N)»(((DTR»$FR(N,1)+1.)»$FR(N,2) 

+.5»(1.-$FR(N,2)))»$FR(N,3)+.25»(1.-$FR(N,3))) 

Where  DTR  = 2 


The  formula  for  the  advance  of  the  variables,  right 
hand  sides  and  relative  errors  are 

$FR(N,1)  = 1.-$FR(N,3)  I 

$V(N,L,1)  = $V(N,L,1)»$FR(N,1)+$V(N,L,3)*$FR(N,3)  I 

RCTURN  ( 

END  , I 

l»)  PRINTOUT 

An  example  of  the  diagnostic  printout  which  may  be 
desired  by  the  user  is  given. This  is  acchieved  by  setting  the 

values  of  LPRI  in  the  calling  sequence  of  AR$MAI  i 

cor  respond ing ly . 

The  first  call  will  result  in  the  following: 

ARIS  EXECUTION  HAS  BEEN  AT  03/03/77  10:55 
PROGRAM  EXECUTION  AT  03/03/77  10:58 

> Table  1 ) 

DEST 
9.811E  05 
3.263E  05 
6.558E  05 
9.504E-09 
9.597E-08 

Table  2) 


VALUE 

PROD 

N-f 

9.902E 

07 

O.OOOE  00 

N0+ 

2.274E 

04 

3.237E  05 

02+ 

7.927E 

04 

6.574E  05 

TE 

9.969E-02 

-4.429E-04 

TI 

1.007E 

00 

6.746E-03 

» - 


I 
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RCT=1 

N-f 

6.574E 

05 

6.264E- 

-10 

FUN 

= 02+ 

RCT=2 

N-f 

3.237E 

05 

3.048E- 

-10 

FUN 

= N0+ 

RCT=3 

N0+ 

3.236E 

05 

1.423E 

01 

FUN 

O.OOOE 

00 

9.493E- 

■13 

fUN 

RCT=4 

02+ 

6.558E 

05 

8.272E 

00 

FUN 

Note  that  only  the  names  appear  which  are  variables  and 
are  involved  in  this  reaction 

DEST.  SOURCE  TERMS 

TE  9.504E-09  FUN  TI  9.597E-08  FUN 
PROD.  SOURCE  TERMS 

TE  9.597E-08  FUN  TI  9.504E-09  FUN 
table  3) 

N4- 


EXT 

9.902E  07 

INT 

9.902E  07 

PROD 

DEST 

TOT 

9.81  IE  05 

1 

6.574E  05 

2 

3.237E  05 

N0+ 

EXT 

2.274E  05 

INT 

2.274E  05 

PROD 

DEST 

TOT 

3.237E  05 

TOT 

3.236E  05 

2 

3.237E  05 

3 

3.236E  05 

3 

O.OOOE  00 

02+ 

EXT 

7.927E  04 

INT 

7.927E  04 

PROD 

DEST 

TOT 

6.574E  08 

TOT 

6.558E  05 

1 

6.574E  05 

4 

6.558E  05 

TE 

EXT 

9.969E-02 

INT 

9.968E-02 

PROD 

DEST 

TOT 

-4.429E-04 

TOT 

9.504E-09 

S« 

9.597E-08 

S« 

9.504E-09 

TI 

EXT 

1.007E  00 

INT 

1.007E  00 

PROD 

DEST 

TOT 

6.743E-03 

TOT 

9.597E-08 

S» 

9.504E-09 

S» 

9.579E-08 

j 

I 

i ^ 


1 
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VII.  END  REMARKS 

The  developed  system  will  make  it  easier  to  include  chemical  reactions  in  large  (hydro) 
codes.  Its  application  will  result  in  a better  understanding  of  the  reactions  (and  their 
coefficients)  involved. 

The  resulting  FORTRAN  program  is  efficient  (vectorized)  and  last  but  not  least  always 
correct  if  the  reactions  are  allright.  The  printed  out  differential  equations  are  integrated.' 
Therefore  there  is  no  need  to  check  this  part  (which  is  the  most  cumbersome)  of  the  code. 

The  ARIS  system  has  been  developed  with  the  ASC  of  NRL  in  mind  but  the  aris  program 
itself  is  easily  transferable  to  other  computers  as  it  is  written  in  standard  FORTRAN,  (some  bit 
handling  parts  may  have  to  be  changed.)  The  ARIS  system  has  been  already  successfully  incor- 
porated in  several  hydrocodes  of  NRL.  It  has  been  compared  with  older  chemical  routines 
(CHEMEQ).  It  gives  the  same  results  in  approximate  the  same  number  of  operations.  For 
very  stiff  equations  the  integrationscheme  (OPT-1)  seems  to  be  superior  to  the  old  CHEMEQ 
(essentially  OPT  — 3). 
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